setwd("/home1/30/jc227089/SRE_inv/ran_outs")

nfiles<-length(grep("SRE_inv_ParSpace", list.files()))
d<-c()
for (i in 1:nfiles){
	fname<-paste("SRE_inv_ParSpace", i, ".RData", sep="")
	load(fname)
	d<-rbind(d, sim.result)
	rm(sim.result)	
}
save(d, file="ParSpaceConcat.RData")
setwd("/home1/30/jc227089/SRE_inv/Figures")
diff<-d[,3]-d[,2]
diff<-ifelse(is.na(diff), 3000, diff)  #treat NAs as 3000 for calculation of mean
tv<-matrix(1:(length(d[,1])/2), nrow=5)
tv<-as.vector(rbind(tv, tv))
d<-cbind(d, diff, tv)
rm(diff, tv)
d2<-cbind(d[which(d[,"spread"]==30),], tapply(d[,"diff"], d[,"tv"], diff))
d3<-aggregate(d2, list(K=d2[,"K"], lambda=d2[,"lambda"], d.cost=d2[,"d.cost"]), mean)


dk40<-subset(d3, d3[,"K"]==40)
dk100<-subset(d3, d3[,"K"]==100)

Y<-c(0, 0.1, 0.5, 0.9)
X<-c(2, 5, 10, 15, 20)
mat40<-matrix(dk40[,"V10"], nrow=length(X))
mat100<-matrix(dk100[,"V10"], nrow=length(X))

png("Par_Space.png", width=12, height=24, units="cm", res=300)
par(mfrow=c(2,1))
contour(X, Y, mat40,ylab="Dispersal cost", xlab="", main=expression(italic("K")*" = 40"))
mtext(text="(a)", side=3, line=1, adj=0)
contour(X, Y, mat100,ylab="Dispersal cost", xlab=expression("Maximum offspring "*italic("R")["max"]), main=expression(italic("K")*" = 100"))
mtext(text="(b)", side=3, line=1, adj=0)
dev.off()


